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We  have  examined  the  Maxwell-Garnett,  inverted  Maxwell-Garnett,  and  Bruggeman  rules  for  evaluation 
of  the  mean  permittivity  involving  partially  empty  cells  at  particle  surface  in  conjunction  with  the 
finite-difference  time-domain  (FDTD)  computation.  Sensitivity  studies  show  that  the  inverted  Maxwell- 
Garnett  rule  is  the  most  effective  in  reducing  the  staircasing  effect.  The  discontinuity  of  permittivity  at 
the  interface  of  free  space  and  the  particle  medium  can  be  minimized  by  use  of  an  effective  permittivity 
at  the  cell  edges  determined  by  the  average  of  the  permittivity  values  associated  with  adjacent  cells.  The 
efficiency  of  the  FDTD  computational  program  is  further  improved  by  use  of  a  perfectly  matched  layer 
absorbing  boundary  condition  and  the  appropriate  coding  technique.  The  accuracy  of  the  FDTD  method 
is  assessed  on  the  basis  of  a  comparison  of  the  FDTD  and  the  Mie  calculations  for  ice  spheres.  This 
program  is  then  applied  to  light  scattering  by  convex  and  concave  aerosol  particles.  Comparisons  of  the 
scattering  phase  function  for  these  types  of  aerosol  with  those  for  spheres  and  spheroids  show  substantial 
differences  in  backscattering  directions.  Finally,  we  illustrate  that  the  FDTD  method  is  robust  and 
flexible  in  computing  the  scattering  properties  of  particles  with  complex  morphological  configurations. 
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1.  Introduction 

The  finite-difference-time  domain  (FDTD)  technique 
is  one  of  the  most  promising  methods  for  solution  of 
the  scattering  properties  of  nonspherical  and  inho¬ 
mogeneous  particles.1-3  The  fundamental  principle 
of  this  method  has  remained  essentially  unchanged 
since  its  inception  more  than  30  years  ago,1  that  is, 
the  electromagnetic  field  is  solved  on  a  discrete  lattice 
on  the  basis  of  time-dependent  Maxwell  equations. 
In  the  FDTD  computation  the  spatial  domain  must  be 
discretized  by  use  of  a  grid  mesh.  It  is  well  known 
that  a  Cartesian  grid  is  more  practical  and  flexible 
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than  a  curvilinear  or  target-conforming  grid  mesh 
because  the  latter  type  of  grid  mesh  is  usually  appli¬ 
cable  only  to  some  specific  geometries.  In  addition,  a 
non-Cartesian  difference  scheme  may  substantially 
complicate  the  criterion  of  computational  stability, 
boundary  conditions,  and  the  specification  of  the  ra¬ 
diation  source.  When  a  particle  with  a  curved  sur¬ 
face  is  defined  over  a  Cartesian  grid,  its  geometry  is 
approximately  represented  by  a  number  of  cubic  cells 
in  a  staircasing  manner.  The  staircasing  effect  can 
potentially  be  large  if  the  dielectric  constants  of  the 
particle  differ  substantially  from  those  of  the  host 
medium.  To  reduce  this  effect  in  the  FDTD  compu¬ 
tation,  Yang  and  Liou2’3  applied  the  Maxwell-Garnett 
rule4  to  evaluate  the  mean  permittivity  of  the  par¬ 
tially  empty  cells  at  the  particle  surface  where  sub¬ 
grid  variations  of  the  optical  constants  exist. 
However,  various  methods  can  be  used  to  compute 
the  mean  permittivity,  and  their  relative  perfor¬ 
mance  requires  further  study. 

In  addition  to  the  staircasing  effect,  a  discontinuity 
of  permittivity  exists  at  a  cell  edge  located  at  the 
interface  of  a  particle  medium  and  free  space,  even  if 
the  grid  mesh  is  target  conforming  (e.g.,  the  Carte¬ 
sian  grid  for  a  cubic  geometry),  a  fact  that  had  not 
been  addressed  in  the  FDTD  implementation.  In 
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the  time-marching  iterations  for  computing  the  near 
field,  permittivity  values  at  cell  edges  are  required. 
Thus  an  appropriate  averaging  procedure  must  be 
used  to  overcome  the  discontinuity  of  permittivity  at 
a  cell  edge  if  the  optical  constants  for  the  four  adja¬ 
cent  cells  are  different.  In  this  study  we  have  im¬ 
proved  the  Cartesian  FDTD  scheme  by  minimizing 
the  staircasing  effect  and  the  discontinuity  of  permit¬ 
tivity  at  the  interface  of  the  particle  medium  and  free 
space.  Because  the  performance  of  the  absorbing 
boundary  is  critical  to  the  accuracy  of  the  FDTD  so¬ 
lution,  we  also  applied  the  newly  developed  perfectly 
matched  layer  (PML)  absorbing  boundary  condi¬ 
tion.5"6  In  Section  2  we  present  three  approaches  to 
evaluating  the  mean  dielectric  constants  for  a  given 
grid  cell  and  suggest  an  appropriate  way  to  compute 
effective  dielectric  constants  at  cell  edges.  To  apply 
the  PML  technique  we  developed  an  equivalent  ver¬ 
sion  of  Berenger’s  PML  equations.6  In  Section  2  the 
program-coding  technique  is  discussed  in  terms  of 
reduction  of  CPU  time  and  memory  requirements. 
Section  3  describes  the  numerical  tests  that  we  car¬ 
ried  out  to  investigate  the  relative  performance  of 
various  approaches  to  determining  the  mean  and  ef¬ 
fective  dielectric  constants  for  grid  cells  and  cell 
edges.  Applications  to  light  scattering  by  non- 
spherical  aerosols  are  described  in  Section  4.  Con¬ 
clusions  are  given  in  Section  5. 

2.  Cartesian  Finite-Difference  Time-Domain  Scheme 


—  Maxwell-Garnett  •  -  -  Inverted  Maxwell-  Garnett  - Bruggeman 


Fig.  1.  Mean  permittivity  evaluated  from  Maxwell-Garnett,  in¬ 
verted  Maxwell-Garnett,  and  Bruggeman  rules  for  grid  cells  com¬ 
posed  of  free  space  and  ice  at  0.55  and  15  pm.  The  refractive 
indices  for  the  two  wavelengths  are,  respectively,  1.311  +  i3. 11  X 
10~9  and  1.571  +  iO.1756. 


cal.  If  we  transpose  the  roles  of  the  two  media  in  the 
Maxwell-Garnett  rule,  we  obtain  the  inverted 
Maxwell-Garnett  rule  as  follows: 

em  +  2e  +  2(1  -  f){em  +  2e)l 

8  =  8  -  .  (2) 

em  +  2e  -  (1  -  f)(em  +  2e) 


A.  Mean  Permittivity  for  Partially  Empty  Cells 
When  a  nonrectangular  particle  is  defined  in  a  Car¬ 
tesian  grid  mesh,  grid  cells  located  near  the  particle 
surface  may  be  partially  empty.  Because  subgrid 
features  cannot  be  specified  in  the  Cartesian  discrete 
lattice,  a  brute-force  approach  is  usually  used  for  de¬ 
termining  the  dielectric  constants  for  these  surface 
cells;  that  is,  a  cell  is  considered  to  be  empty  (with 
unity  permittivity)  if  half  of  the  cell  volume  is  outside 
the  particle.  Otherwise,  the  dielectric  properties  of 
the  particle  are  assigned  to  the  cell.  This  approach 
produces  a  sharp  truncation  in  terms  of  dielectric 
properties  in  representing  the  scattering  particle  by  a 
step-by-step  approximation,  and  the  pertinent  stair¬ 
casing  effect  may  be  potentially  large.  Thus  it  is 
necessary  to  evaluate  the  mean  dielectric  constant  for 
the  partially  empty  cells. 

The  Maxwell-Garnett4  and  Bruggeman7  rules  have 
commonly  been  applied  to  the  evaluation  of  the  mean 
dielectric  constant  associated  with  inhomogeneity. 
For  a  given  volume  containing  two  dielectric  media 
(host  and  inclusion),  the  mean  dielectric  constant 
given  by  the  Maxwell-Garnett  rule  is 

__  [e  +  2em  +  2/’(e  +  2ejl 

^  Ln  if)  />/  ,  f>  \  ’  (1) 

[  e  +  2em  -  f(e  +  2eJ  J 

where  em  and  e  are  the  permittivities  of  the  host  and 
the  inclusion  media,  respectively,  and/7 is  the  volume 
fraction  of  the  inclusion.  It  can  be  seen  that  the 
roles  of  host  and  inclusion  in  Eq.  (1)  are  not  recipro¬ 


The  Bruggeman  rule  can  be  expressed  in  the  follow¬ 
ing  general  form: 


1 

Av 


e(r)  -  e 
e(r)  +  2e 


do  =  0, 


(3) 


where  Av  is  the  spatial  region  for  accounting  for  the 
inhomogeneity  of  dielectric  properties.  If  only  two 
media  are  considered,  Eq.  (3)  reduces  to 


_  8  —  8  Em  8 

f— Oz  +  a-n^v^  =  o- 

e  +  2e  em  +  2e 


(4) 


Note  that  the  roles  of  the  host  and  inclusion  media 
are  reciprocal  in  the  Bruggeman  equation.  Equa¬ 
tion  (4)  is  nonlinear  and  complex,  and  its  solution  is 
not  straightforward.  For  practical  calculations,  we 
have  derived  an  efficient  iterative  scheme  for  the  so¬ 
lution  of  Eq.  (4),  given  by 


_  =  eem  +  2  [/e  +  (1  -/)£„]£„ 

e"+1  /£m  +  (1  -/■)£  +  2e„ 

With  an  initial  value  of  £0  =  fe  +  (1  —  f)em,  a  conver¬ 
gent  solution  can  be  obtained  with  iterations  in  fewer 
than  ten  steps. 

Figure  1  shows  the  mean  permittivity  values  eval¬ 
uated  from  Eqs.  (1),  (2),  and  (4)  for  a  cell  composed  of 
free  space  and  ice  at  0.55-  and  15-p.m  wavelengths. 
The  refractive  indices  of  ice  for  these  two  wave¬ 
lengths  are  1.311  +  (3.11  X  10~9  and  1.571  + 
(0.1756,  respectively.  Note  that  the  complex  per- 
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mittivity  (er  +  iz,)  is  related  to  the  complex  refractive 
index  ( m +  im,)  as  follows: 


er  =  mr2  —  m2,  zt  =  2mrmi.  (6) 


From  Fig.  1  it  is  evident  that  the  differences  in  per¬ 
mittivity  evaluated  from  the  three  schemes  are  very 
small  for  the  optically  thin  case  at  0.55  pm.  Because 
the  absorption  of  ice  at  0.55  pm  is  extremely  small,  it 
does  not  make  a  noticeable  difference  in  the  evaluation 
of  the  mean  imaginary  part  of  the  permittivity.  For 
the  15-pm  wavelength,  however,  the  imaginary  part  of 
the  refractive  index  is  quite  large.  In  this  case,  the 
absorption  properties  computed  from  the  FDTD  tech¬ 
nique  depend  on  the  approach  used  to  evaluate  the 
mean  permittivity.  The  mean  permittivity  evaluated 
by  the  Bruggeman  rule  lies  between  those  evaluated 
from  the  Maxwell-Garnett  and  inverted  Maxwell- 
Garnett  rules.  Results  computed  from  Bruggeman’s 
equation  converge  to  those  determined  from  Maxwell- 
Garnett’s  equation  if  the  cells  are  nearly  empty,  and 
they  converge  to  the  values  evaluated  from  the  in¬ 
verted  Maxwell-Garnett  rule  when  ice  dominates  the 
cells.  However,  the  overall  patterns  of  variation  of 
mean  permittivity  versus  fraction  of  host  medium  are 
similar  for  the  three  approaches. 


B.  Effective  Permittivity  at  Cell  Edges 
The  finite-difference  analog  of  Maxwell’s  equations 
has  a  staggered  form;  that  is,  the  electric  field  is 
evaluated  at  the  cell  edges  at  time  step  n,  whereas 
the  magnetic  field  is  evaluated  at  the  center  of  cell 
faces  at  time  step  n  +  1/2.  Thus,  for  time-marching 
iterative  computation  of  the  electric  field,  the  permit¬ 
tivity  values  at  the  cell  edges  must  be  known.  It  is 
straightforward  to  specify  the  permittivity  value  at  a 
cell  edge  if  the  four  adjacent  cells  have  the  same 
dielectric  properties.  However,  there  is  a  disconti¬ 
nuity  or  singularity  of  permittivity  at  the  cell  edge  if 
the  four  adjacent  cells  are  heterogeneous.  In  this 
case  it  is  neither  mathematically  suitable  nor  physi¬ 
cally  correct  to  assume  indiscriminately  that  the 
value  of  the  permittivity  at  the  cell  edge  is  the  value 
associated  with  one  of  the  adjacent  cells,  although 
this  approach  has  commonly  been  employed  for  nu¬ 
merical  realization  of  the  FDTD  technique.  To  de¬ 
rive  an  appropriate  value  of  permittivity  at  the  cell 
edge  we  begin  with  the  discretization  of  Maxwell’s 
equations.  Consider  the  z  component  of  the  E  field 
as  an  example.  The  equation  that  governs  the  tem¬ 
poral  variation  of  an  electric  field  is  given  by 


er(r)  dE(r,  t) 

C  dt 


+  kz i(r)E(r,  t ) 


V  X  H(r,  t ), 


(7) 


Fig.  2.  Geometry  for  a  cell  edge  and  four  adjacent  cells  in  a 
Cartesian  grid. 


sion  in  Eq.  (7)  over  the  domain  enclosed  by  the  dotted 
lines  shown  in  Fig.  2  leads  to 

(M-l)Ax  f(J+l)Ay 

[V  X  H(r,  Gldxdy 

/Ax  J  J\y 

=  A x{Hx[r(I  +  1/2,  J,  K),  t]  -  HX[?(I  +  1/2, 

J+1,K),  t]}  +  A y{Hy[r{I  +  1  ,J+  1/2,  K),  t ] 

-  Hy[r(I,  J  +  1/2,  K),  t]\,  (8) 

where  r{I  +  1/2,  J,  K)  represents  the  position  vector 
for  point  [(/  +  l/2)Ax,  JAy,  kA z\  It  should  be 
pointed  out  that  the  discretization  in  Eq.  (8)  does  not 
produce  any  truncation  error.  Furthermore,  inte¬ 
grating  the  left-hand  side  of  Eq.  (7)  yields 


(7+l)Ax  /'(=/+!)  Ay 


er(r)  dEz(r,  t ) 


dt 


+  kZi{r)Ez(r ,  t) 

.  C  cil 

I  Ax  *JAy 

z r(I  +  1/2,  J  +  1/2,  K )  d 


dxd y 


dt 


*(/+l) Ax  /* 

JJ 

1 

IAx  v 

J+  1/2,  K) 

E r(I  +  1/2,  J  +  1/2,  K)  d 


dt 


+  kz  i(I  +  1/2, 
Ez(r,  t)dxdy 
+  kz,i(I  +  1/2, 


J  +  1/2,  K) 


Ez[r{I  +  1/2,  J  +  1/2,  K),  t\AxAy,  (9) 


where  £,.(/  +  1/2,  J  +  1/2,  K)  and  e t{I  +  1/2,  J  +  1/2, 
K)  are  the  effective  permittivities  evaluated  at  the 
cell  edge.  From  Eq.  (8)  and  relation  (9)  we  obtain 
the  following  difference-differential  equation: 

er[7(/  +  1/2,  J  +  1/2,  K)] 


dEz[f(I  +  1/2,  J  +  1/2,  K),  t ] 


dt 


where  c  is  the  speed  of  light  in  vacuum,  k  =  w/c,  where 
to  is  the  angular  frequency  of  the  electromagnetic 
wave,  and  the  real  and  imaginary  parts  of  permittivity 
can  be  computed  from  Eqs.  (6).  Let  the  cell  center  be 
denoted  (x,  y,  z)  =  (/Ax,  JAy,  KAz ),  where  Ax,  Ay,  and 
Az  are  the  cell  dimensions  along  the  three  coordinate 
axes.  Integration  of  the  z  component  of  the  expres- 


+  kz ,[r(/  +  1/2,  J  +  1/2,  Kj\Ez[HI  +  1/2, 

J  +  1/2,  K),  t ] 

=  {Him  +  1/2,  j,  K),  t]  -  Hx\r(I  +  1/2, 

J  +  1,  K),  t]}/Ay  +  {Him  +  1,  j  +  1/2,  K ),  t] 

-  Him,  j  +  1/2,  K),  tl/ Ax.  (10) 
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The  temporal  derivative  in  Eq.  (10)  can  be  discretized 
in  a  straightforward  manner  by  use  of  either  a  leap¬ 
frog  or  an  exponential  scheme.  Note  that  the  dis¬ 
cretization  in  time  will  not  affect  the  effective 
permittivity  at  the  cell  edge  because  it  is  a  time- 
independent  variable.  To  evaluate  the  effective  per¬ 
mittivity  we  let  the  complex  permittivity  values  for 
the  four  adjacent  cells  be  e1  =  Erl  +  is, , ,  e2  =  er2  +  ie(2, 
£3  =  e„3  +  iei3,  and  e4  =  er4  +  ie(>1.  It  follows  that 


sions  in  Maxwell’s  equations.  For  a  systematic  pre¬ 
sentation  of  this  study,  here  we  recapitulate  the  PML 
boundary  equations  expressed  in  a  format  that  con¬ 
tains  diffusion  terms. 

The  diffusions  are  imposed  on  the  electromagnetic 
field  components  that  are  parallel  to  the  boundary. 
In  practice,  the  field  must  to  be  separated  into  two 
components,  one  parallel  and  one  perpendicular  to 
the  boundary,  that  is, 


£r(/  +  1/2,  J  +  1/2,  K)  — 


+  £/-4 


Li 

%  (1+1/2)  Ax 

'(J+l/2)Ay 

+  L-2 

»(/+  l)Ax 

'(</+ 1/2)  Ay 

+  L-3 

*(/+l)Ax 

i 

/A*  ‘ 

JAy  * 

(1+1/2) Ax  ‘ 

JAy  * 

(1+1/2) Ax  J 

C  (1+1/2)  Ax 

•(J+l)Ay 

lAx  ‘ 

(J+l/2)Ay_ 

Ez(r,  Dd.vdv 


Ez(r,  Ddxdy  «  (eri  +  er2  +  er3  +  Erd/4.  (11) 


“J Ay 


The  error  for  the  approximation  in  relation  (11)  is 
small  because  the  tangential  component  of  the  elec¬ 
tric  field  is  continuous  across  the  interface,  as  re¬ 
quired  by  the  electromagnetic  boundary  condition. 
Similarly,  we  have 


(E„  Ey,  Ez)  =  [(Exy  +  EJ,  (Eyx  +  Eyz),  (Ezx  +  Ezy)\ 

(13a) 

(H„  Hy,  Hz)  =  [(Hxy  +  HJ,  (Hyx  +  HyJ,  (Hzx  +  Hzy)], 

(13b) 


£i(I+  1/2,  J+  1/2,  K)  «  (ea  +  ei2  +  e,:i  +  e,4)/4.  (12)  The  PML  equations  for  the  Ex  and  Hx  components,  for 

example,  are 


The  preceding  procedure  of  averaging  the  cell  per¬ 
mittivity  is  mandatory  for  any  particle  geometry  de¬ 
fined  in  a  Cartesian  grid,  even  if  the  staircasing 
approximation  for  defining  the  particle  geometry  is 
absent.  For  example,  a  cubic  particle  can  be  well 
defined  in  the  Cartesian  grid  without  the  staircasing 
approximation.  However,  the  cell  edges  at  the  par¬ 
ticle  surface  are  surrounded  by  both  empty  cells  and 
dielectric  medium  cells,  a  condition  for  which,  to 
avoid  discontinuity,  an  average  of  the  permittivity  is 
required. 

C.  Perfectly  Matched  Layer  Absorbing  Boundary 
Condition 

The  performance  of  the  absorbing  boundary  condition 
is  critical  to  the  accuracy  of  the  FDTD  technique. 
Various  analytical  absorbing  boundary  equations 
have  been  developed  in  the  past  two  decades  (see  Ref. 
8  and  the  papers  cited  therein).  These  boundary 
equations  usually  require  a  substantial  white  space 
between  the  boundary  and  the  scattering  target, 
thereby  reducing  the  efficiency  of  the  FDTD  method 
in  terms  of  CPU  and  memory  requirements.  Be- 
renger’s  PML5-6  has  been  widely  applied  to  various 
FDTD  implementations.  Numerical  experiments 
have  shown  that  the  spurious  reflection  associated 
with  the  PML  boundary  condition  is  approximately  3 
orders  of  magnitude  less  than  that  produced  by  an 
analytical  boundary  equation.9  The  PML  technique 
has  been  successfully  applied  to  light  scattering  by 
ice  crystals.10  Yang  and  Liou11  have  shown  that  the 
absorption  of  the  PML  boundary  condition  is  funda¬ 
mentally  due  to  the  introduction  of  artificial  diffu- 


exp[— Ty(y)t]  d  d(Hzx  +  Hzy) 

{exp[Ty(y)t]F,y}  = 


dt 


exp[-Tz(z)f]  d 
C  dt 


{exp[T z(z)t]Exz\  =  - 


dy 

(14a) 

d(Hyx  +  Hyz) 


dz 


(14b) 


exp[-Ty(y)t]  a  d(Ezx  +Ezy) 

{exp[Ty(y)f]Hyy}  =  - 


dt 


dy 


(14c) 


exp[-Tz(z)f]  d 
C  dt 


1  r  7  uicr  \  d(Eyx  +  EyJ 

{exp[T z(z)t]HJ  = - - - — , 

dz 


(14d) 


where  t y(y)  and  ~Jz)  are  zero,  except  in  the  boundary 
layers  perpendicular  to  the  y  and  z  axes.  The 
present  form  for  the  PML  boundary  equations  is 
equivalent  to  that  given  by  Berenger.6  In  practical 
computations,  the  diffusion  parameters  t  (y)  and 
t Z{z)  can  be  specified  from  zero  at  the  interface  of  free 
space  and  the  artificial  diffusion  medium  to  their 
maximum  values  at  the  outermost  layer.  For  exam¬ 
ple,  t y(y)  can  be  defined  as 

T y(y)  =  rTy,max[(y  -y0)/DJ,  (15) 

where  (y  -  y0)  is  the  distance  of  a  grid  point  from  the 
interface  of  free  space  and  the  PML  medium,  D  = 
LAy  is  the  thickness  of  the  artificial  diffusion  medium 
for  the  boundary  perpendicular  to  the  y  axis,  and p  is 
usually  selected  to  be  in  the  range  2-3.  The  param- 
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eter  jy  max  can  be  specified  in  terms  of  the  reflectance 
of  the  boundary  for  normal  incidence  as  follows: 

A, max  =  ln[^(0°)]C-  (16) 

where  R(0°)  is  the  boundary  reflection  factor  for  nor¬ 
mal  incidence.  The  mean  absorption  must  be  taken 
into  account  for  each  cell  distance  in  discrete  compu¬ 
tations.  Thus  the  following  two  mean  values  for  the 
electric  and  the  magnetic  fields  can  be  used: 


T y(J) 


(J+  l/2)Ax 


Ty(y)dy 


J-1/2A.X 


W  (J  +  l/2f+1  -  (J  -  l/2f+1 

p  +  1  Lp+1 


E  field, 
(17a) 


1  P(J+  l)Ay 

t y(J  +  1/2)  =  —  Ty(y)dy 

y  JjAy 

_  w  (j+ir1-^1 

p  +  1  Lp+1 


H  field. 

(17b) 


In  this  study  we  use  the  PML  absorbing  boundary 
condition  with  seven  layers  of  artificial  diffusion  me¬ 
dium  that  are  separated  from  the  scattering  target  by 
eight  layers  of  white  space.  According  to  the  sensi¬ 
tivity  study  of  Lazzi  and  Gandi,12  we  select  p  =  2.5 
and  R{ 0°)  =  5  X  1(T6 


D.  Program  Coding  Improvement 
Advances  in  the  FDTD  algorithms  coupled  with  the 
increasing  power  of  computers  now  make  it  possible 
to  apply  this  method  to  size  parameters  up  to  20  on  an 
ordinary  workstation.  For  particles  with  refractive 
index  close  to  unity,  FDTD  computation  has  been 
carried  out  on  a  supercomputer  for  size  parameters 
up  to  40. 10  The  grid  size  can  be  much  larger  in  an 
optically  thin  case  than  in  an  optically  thick  case. 
The  huge  CPU  and  memory  requirements  of  the 
FDTD  technique  are  still  a  major  concern  in  practical 
applications.  However,  a  substantial  reduction  in 
computer  resources  can  be  achieved  by  proper  pro¬ 
gram  coding. 

During  the  course  of  computational  experimenta¬ 
tion  we  found  that  the  efficiency  of  the  FDTD  tech¬ 
nique  can  be  greatly  improved  if  the  computational 
program  is  coded  in  fortran90  rather  than  in  for- 
tran77.  In  particular,  the  array-oriented  feature  in 
fortran90  can  substantially  reduce  the  CPU  time  for 
initialization  and  updating  of  the  huge  arrays  in¬ 
volved  in  the  time-marching  iterations  of  the  near 
field  and  in  the  Fourier-transform  calculations  that 
yield  the  field  values  in  the  frequency  domain.  In 
addition,  fortran90  permits  dynamic  storage  alloca¬ 
tion.  With  this  feature  the  temporal  working  arrays 
that  specify  particle  constants  and  the  PML  coeffi¬ 
cients  can  be  deallocated  immediately  as  long  as  they 


are  not  required  for  further  computation.  The  effi¬ 
ciency  of  the  FDTD  computational  program  can  also 
be  improved  by  use  of  various  intrinsic  functions  that 
are  permitted  in  fortran90.  Yang  and  Liou3  have 
formulated  an  algorithm  to  map  the  near  field  (in  the 
frequency  domain)  to  the  corresponding  far  field  on 
the  basis  of  a  volume  integral  of  the  electric  field 
inside  the  scattering  particle.  With  the  dynamic 
memory  allocation  feature  of  fortran90,  three  one¬ 
dimensional  arrays  can  be  allocated  with  the  exact 
array  sizes  in  the  computations  to  handle  the  three 
Cartesian  components  of  the  E  field  with  only  non¬ 
empty  cells.  Moreover,  the  FDTD  computation  can 
be  largely  vectorized  with  these  one-dimensional  ar¬ 
rays  in  the  Fourier-transform  procedure  for  the  field 
values  in  the  frequency  domain  and  in  mapping  the 
near  file  to  the  far  field. 

3.  Numerical  Tests  and  Discussions 

We  have  improved  the  FDTD  program  for  light  scat¬ 
tering  by  dielectric  particles  developed  by  Yang  and 
Liou3  by  updating  the  absorbing  boundary  condition 
with  the  PML  technique  as  well  as  by  coding  the 
program  in  fortran90.  For  a  size  parameter  of  10, 
the  CPU  requirement  is  reduced  by  a  factor  of  2  and 
the  memory  requirement  is  decreased  by  approxi¬ 
mately  10%.  The  improvement  is  more  substantial 
for  larger-sized  parameters.  With  the  improved  pro¬ 
gram,  we  conducted  sensitivity  computations  in  con¬ 
junction  with  the  various  schemes  for  evaluating 
mean  and  effective  permittivities.  Because  the 
FDTD  technique  does  not  give  preferential  treatment 
to  any  specific  particle  shape  (with  the  exception  of 
rectangular  targets  in  a  Cartesian  grid),  a  canonical 
study  based  on  spheres  can  be  a  representative  test  of 
the  technique’s  performance.  For  this  reason,  the 
present  computations  focus  on  the  scattering  proper¬ 
ties  of  ice  spheres  at  the  15- qm  infrared  wavelength 
at  which  the  refractive  index  differs  substantially 
from  unity.  A  spherical  particle  defined  in  a  Carte¬ 
sian  grid  lattice  is  a  pseudosphere  with  step-by-step 
surface  roughness.  The  magnitude  of  the  surface 
roughness  or  staircasing  effect  is  proportional  to  the 
grid  cell  size  and  inversely  proportional  to  the  parti¬ 
cle  dimension.  It  follows  that  the  present  sensitivity 
study  focuses  on  two  representative  cases  with  size 
parameters  of  1  and  10. 

Figure  3  shows  a  comparison  of  the  phase  functions 
of  ice  spheres  computed  from  Mie  theory  and  from  the 
FDTD  technique  for  a  size  parameter  of  1.  In  the 
latter  technique  we  applied  the  Maxwell-Garnett,  in¬ 
verted  Maxwell-Garnett,  and  Bruggeman  rules  to 
evaluate  the  mean  permittivity  for  the  cells  near  the 
particle  surface.  In  addition,  we  used  an  average  of 
the  permittivity  given  by  relations  (11)  and  (12)  to 
define  the  effective  permittivity  at  cell  edges  if  the 
adjacent  cells  are  different  in  dielectric  constants. 
Three  grid  sizes,  As  =  X/20,  X/25,  X/30,  were  se¬ 
lected,  where  X  is  the  wavelength  in  vacuum  and  As 
(=Ax  =  A y  =  A z)  is  the  grid  cell  size.  From  Fig.  3  it 
can  be  seen  that  the  inverted  Maxwell-Garnett  rule 
performs  best  of  the  three  approaches  to  evaluating 
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Scattering  Angle  (deg.) 

Fig.  3.  Comparison  of  FDTD  and  Mie  results  for  the  phase  func¬ 
tions  of  ice  spheres  at  15  p.m  for  a  size  parameter  of  1.  Three 
approaches  are  used  for  evaluating  the  mean  permittivity  for  par¬ 
tially  empty  cells  located  near  the  particle  surface. 


mean  permittivity.  In  addition,  it  is  evident  that  the 
accuracy  of  the  FDTD  results  is  sensitive  to  the  ratio 
of  grid  size  to  incident  wavelength.  The  effect  of  grid 
size  on  the  accuracy  of  the  FDTD  solution  has  been 
discussed  by  Yang  and  Liou.3 

Figure  4  is  similar  to  Fig.  3,  except  that  the  size 
parameter  is  now  10.  For  this  size  parameter  the 
effect  of  surface  roughness  as  a  result  of  the  staircas¬ 
ing  approximation  of  the  particle  shape  is  10  times 
smaller  than  that  shown  in  Fig.  3  for  size  parameter 
1.  Note  that  the  magnitude  of  the  roughness  effect 
is  proportional  to  A s/r,  where  r  is  the  radius  of  the 
sphere.  Thus  the  differences  in  the  results  obtained 
from  the  three  methods  applied  to  evaluation  of  the 
mean  permittivity  for  the  partially  empty  cells  are 
negligible  for  a  large  size  parameter.  Parameter 
X/As  is  a  key  factor  for  determining  the  accuracy  of 
the  FDTD  technique,  as  is  evident  from  Fig.  3. 

Table  1  lists  the  extinction  efficiencies  and  single¬ 
scattering  albedos  associated  with  the  phase  func¬ 
tions  for  X/As  =  30  shown  in  Figs.  3  and  4. 
Comparisons  of  absolute  and  relative  errors  of  the 
FDTD  solution  with  the  Mie  results  are  listed  in  the 
table.  For  size  parameter  1,  the  inaccuracy  of  the 
Maxwell-Garnett  rule  is  approximately  twice  that  as¬ 
sociated  with  the  inverted  Maxwell-Garnett  rule. 
The  accuracy  of  the  Bruggeman  rule  lies  between 
those  of  the  two  other  rules,  as  is  expected  from  Fig. 


Fig.  4.  Same  as  Fig.  3  but  for  a  size  parameter  of  10. 


1.  For  size  parameter  10,  the  effect  of  staircasing  is 
small  because  the  roughness  of  the  pseudosphere  de¬ 
fined  in  the  Cartesian  grid  is  insignificant  in  compar¬ 
ison  with  the  particle  dimension.  Again,  for  this  size 
parameter,  the  inverted  Maxwell-Garnett  rule  gives 
the  best  results. 

To  study  the  effect  of  the  effective  permittivity  at 
cell  edges,  we  have  designed  three  schemes.  In 
scheme  1  we  apply  the  inverted  Maxwell-Garnett 
rule  to  evaluate  the  mean  permittivity  for  the  par¬ 
tially  empty  cells  located  at  the  particle  surface.  In 
addition,  we  obtain  the  effective  permittivity  at  a  cell 
edge  in  the  time-marching  iteration  of  the  E  field  by 
averaging  the  permittivity  values  associated  with  the 
four  adjacent  cells.  In  scheme  2,  a  cell  is  considered 
to  be  empty  if  50%  of  its  volume  is  outside  the  parti¬ 
cle,  and  a  unity  permittivity  is  assigned  to  the  cell; 
otherwise,  the  permittivity  of  the  particle  is  applied 
to  the  cell.  In  this  scheme,  averaging  of  the  permit¬ 
tivity  is  not  applied  to  the  cell  edge;  instead,  the 
particle  permittivity  is  assigned  to  the  cell  edge  if  any 
one  of  the  adjacent  cells  is  nonempty.  If  all  four 
adjacent  cells  are  empty,  a  unity  permittivity  is  as¬ 
signed  to  the  cell  edge.  The  procedure  for  construct¬ 
ing  the  particle  shape  in  scheme  3  is  the  same  as  for 
scheme  2.  However,  in  scheme  3  the  averaging  pro¬ 
cedure  is  carried  out  for  the  cell  edges  on  the  basis  of 
relations  (11)  and  (12). 

Figure  5  shows  the  phase  functions  computed  with 
the  these  three  schemes  for  size  parameter  1  and 
their  relative  errors  in  comparison  with  the  Mie  re- 
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Table  1.  Extinction  Efficiencies  Qe  and  Single-Scattering  Albedos  to  Corresponding  to  the  Phase  Functions  Shown  in  Figs.  3  and  4  for  K/As  =  30“ 


Qe  (FDTD) 

A  Qe 

A QJQe  (Mie) 

(%) 

*  (FDTD) 

Act) 

Aoi/u)  (Mie) 

(%) 

X  =  1  [Qe  (Mie)  =  0.7342,  * 

(Mie)  =  0.3699] 

Maxwell-Garnett 

0.7931 

0.0589 

8.02 

0.3799 

0.0100 

2.70 

Inverted  Maxwell-Garnett 

0.7683 

0.0341 

4.64 

0.3757 

0.0058 

1.57 

Bruggeman 

0.7822 

0.0480 

6.54 

0.3772 

0.0073 

1.97 

X  =  10  [Qe  (Mie)  =  2.4177, 
w  (Mie)  =  0.4943] 

Maxwell-Garnett 

2.3904 

-0.0273 

-1.13 

0.4987 

0.0044 

0.89 

Inverted  Maxwell-Garnett 

2.3927 

-0.0250 

-1.03 

0.4983 

0.0040 

0.81 

Bruggeman 

2.3915 

-0.0262 

-1.08 

0.4986 

0.0043 

0.87 

“Also  listed  are  the  values  of  absolute  error  and  relative  error  for  the  two  quantities.  Note  that  A Qe  =  Qe  (FDTD)  -  Qe  (Mie)  and  Aoo  = 
to(FDTD)  -  to(Mie). 


suits.  We  selected  X./A s  =  25  for  the  computations. 
It  is  evident  that  scheme  2  produces  the  largest  er¬ 
rors,  particularly  for  scattering  angles  larger  than 
120°.  The  signs  of  the  variations  of  relative  error 
versus  scattering  angle  are  similar  for  schemes  2  and 
3  because  we  used  the  same  procedure  in  construct¬ 
ing  the  particle  shape.  However,  using  the  averag¬ 
ing  procedure  to  determine  the  effective  permittivity 
for  cell  edges  substantially  improved  the  accuracy  of 
the  FDTD  technique,  as  is  shown  by  comparison  of 
the  results  computed  from  the  three  schemes.  We 
note  that  the  relative  error  in  forward  scattering  is 
larger  in  schemes  2  and  3  than  in  scheme  1.  The 
magnitude  of  the  phase  function  in  the  forward¬ 
scattering  direction  is  approximately  two  times 
larger  than  in  backscattering.  Thus,  in  terms  of  to¬ 
tal  scattered  energy,  scheme  1  is  also  accurate. 

Figure  6  is  the  same  as  Fig.  5,  except  that  the  size 
parameter  is  10.  In  this  case  the  staircasing  approx¬ 
imation  is  small  in  comparison  with  the  particle  di¬ 
mension.  Thus  the  method  of  evaluating  the  mean 


0  60  120  180  0  60  120  180  0  60  120  180 
Scattering  Angle  (deg.) 

Fig.  5.  Comparison  of  the  performance  of  three  schemes  to  ac¬ 
count  for  the  discontinuity  of  permittivity  at  the  interface  of  free 
space  and  the  particle  medium  in  conjunction  with  the  FDTD 
computation  of  phase  function  for  an  ice  sphere  with  a  size  param¬ 
eter  of  1. 


dielectric  constants  for  the  partially  empty  cells  lo¬ 
cated  near  the  particle  surface  becomes  less  impor¬ 
tant.  For  this  reason  the  errors  in  scheme  1  are  of 
the  same  order  as  those  in  scheme  3.  However,  the 
maximum  error  produced  by  scheme  2  is  approxi¬ 
mately  two  times  larger  than  those  produced  by 
schemes  1  and  3.  In  particular,  scheme  2  produces  a 
40%  relative  error  at  the  20°  scattering  angle  where 
a  resonant  minimum  is  shown.  Evidently,  the 
method  with  which  the  effective  dielectric  constant  at 
cell  edges  is  evaluated  can  substantially  affect  the 
performance  of  the  FDTD  technique.  It  affects 
mainly  the  locations  at  the  interface  of  free  space  and 
the  particle  medium,  regardless  of  the  size  parame¬ 
ter.  Physically,  the  accuracy  of  the  electromagnetic 
boundary  condition  (i.e.,  the  tangential  components 
of  the  E  field  and  the  normal  components  of  the  II 
field  are  continuous  at  medium  interface)  involved  in 
the  time-marching  iteration  of  the  near  field  in  the 
time  domain  depends  substantially  on  the  correct¬ 
ness  of  the  dielectric  constants  defined  at  the  inter¬ 
face.  Thus  it  is  critical  to  evaluate  properly  the 
dielectric  constant  at  the  cell  edges. 

Table  2  lists  the  FDTD  solutions  and  the  associated 
errors  for  the  extinction  efficiencies  and  single- 
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Table  2.  Extinction  Efficiencies  Qe  and  Single-Scattering  Albedos  w  Corresponding  to  the  Phase  Functions  Shown  in  Figs.  5  and  6a 


Scheme  Number 

Qe  (FDTD) 

Q? 

<1 

\QJQe  (Mie)  (%) 

*  (FDTD) 

Aw 

Aw/w  (Mie)  (%) 

X  =  1  [Qe  (Mie)  =  0.7342, 

a>  (Mie)  =  0.3699] 

1 

0.7813 

0.0471 

6.42 

0.3814 

0.0115 

3.11 

2 

0.8526 

0.1184 

16.12 

0.4230 

0.0531 

14.35 

3 

X  =  10  \Q„  (Mie)  =  2.4177, 

0.8411 

0.1069 

14.56 

0.3945 

0.0246 

6.65 

<o(Mie)  =  0.4943] 

1 

2.3793 

-0.0384 

-1.59 

0.4994 

0.0051 

1.03 

2 

2.2734 

-0.1443 

-5.97 

0.5022 

0.0079 

1.60 

3 

2.3768 

-0.0407 

-1.68 

0.5003 

0.0060 

1.20 

“Also  listed  are  the  values  of  absolute  error  and  relative  error.  Note  that  A Qe  =  Qe  (FDTD)  —  Qe  (Mie)  and  Aw  =  w  (FDTD)  —  w  (Mie). 


scattering  albedos  that  correspond  to  the  phase  func¬ 
tions  shown  in  Figs.  5  and  6.  It  is  evident  that 
substantial  inaccuracies  can  occur  if  the  effective  per¬ 
mittivity  values  at  the  cell  edges  are  not  properly 
evaluated.  The  inaccuracies  caused  by  inappropri¬ 
ate  treatment  of  the  permittivity  at  the  cell  edges  is 
more  significant  in  the  computation  of  extinction  ef¬ 
ficiency  than  in  the  computation  of  the  single¬ 
scattering  albedo.  This  is  so  because  the  errors  may 
lead  to  an  overestimation  or  underestimation  of  the 
extinction  and  scattering  cross  sections  in  a  consis¬ 
tent  manner;  thus  the  errors  can  be  largely  canceled 
in  the  single-scattering  albedo  calculations.  From 
Figs.  5  and  6  and  Table  2  it  can  be  noted  that  scheme 
1  is  the  most  accurate  of  the  three  schemes.  Thus 
we  shall  use  scheme  1  in  the  rest  of  this  study. 

To  demonstrate  the  accuracy  of  the  present  im¬ 
proved  FDTD  computational  program  we  show  in 
Fig.  7  the  phase  functions  for  an  ice  sphere  at  11-  and 
15-|im  wavelengths.  The  refractive  indices  for  the 
two  wavelengths  are  1.0925  +  i0.248  and  1.571  + 
iO.1756,  respectively.  Because  the  magnitude  of  the 
ice’s  refractive  index  at  the  15-p.m  wavelength  is  ap¬ 
proximately  1.3  times  larger  than  that  at  11  pm,  the 


Scattering  Angle  (deg.)  Scattering  Angle  (deg.) 

Fig.  7.  Comparison  of  FDTD  and  Mie  results  for  the  phase  func¬ 
tion  of  ice  spheres  by  use  of  grid  resolutions  X/As  =  25,  30,  respec¬ 
tively,  for  optically  thin  (1.0925  +  (0.248)  and  thick  (1.5710  + 
iO.1756)  cases. 


resolution  of  X  /A s  =  25  at  the  latter  wavelength  is 
equivalent  to  a  resolution  of  Xp/A s  =  32.5  at  the 
former,  where  M  is  the  wavelength  within  the  parti¬ 
cle  medium.  For  this  reason  we  select  X/A s  =  25, 
30,  respectively,  for  the  11-  and  15-p.m  wavelengths. 
The  relative  errors  of  the  phase  functions  computed 
from  the  FDTD  technique  are  less  than  10%.  The 
error  maxima  are  noted  mainly  at  the  angular  loca¬ 
tions  that  correspond  to  the  resonant  minima  and 
backscattering.  Finally,  it  should  be  pointed  out 
that  the  grid  size  must  be  less  than  1/20  of  the  wave¬ 
length  inside  the  particle  to  prevent  substantial  com¬ 
putational  dispersion  in  the  FDTD  implementation. 

4.  Application  to  Light  Scattering  by  Aerosols 

The  shapes  of  aerosols  are  diverse,  ranging  from 
quasi-spheres  to  highly  irregular  geometries.  Min¬ 
eral,  dustlike,  and  soil  aerosols  are  largely  nonspheri- 
cal  particles  with  mean  aspect  ratios  substantially 
different  from  unity.13-15  Knowledge  of  the  funda¬ 
mental  scattering  and  absorption  properties  of  aero¬ 
sols  is  essential  to  the  development  of  a  physically 
based  and  accurate  retrieval  algorithm  for  aerosol 
optical  depth.16  Mishchenko  et  aZ.17-18  have  demon¬ 
strated  that  the  effect  of  nonsphericity  on  the  optical 
properties  of  aerosols  has  the  potential  to  be  great 
and  that  the  equivalent  spherical  approximation  for 
the  single-scattering  properties  of  nonspherical  par¬ 
ticles  based  on  Mie  theory  is  physically  inappropriate 
and  often  misleading.  Using  a  mixture  of  prolate 
and  oblate  spheroids  for  aerosols,  Mishchenko  et  al.17 
showed  that  the  errors  that  are  due  to  neglect  of 
particle  nonsphericity  are  much  larger  than  those 
that  stem  from  measurement  errors  and  can  easily 
exceed  100%  in  the  retrieval  of  aerosol  optical  depth. 
In  what  follows,  we  investigate  the  deviation  of  the 
phase  functions  computed  for  convex  and  concave 
aerosols  based  on  the  FDTD  technique  from  those 
computed  for  spheres  and  spheroids. 

The  sizes  of  convex  and  concave  aerosols  can  be 
specified  in  terms  of  their  peripheral  spheroids.  To 
construct  a  convex  geometry,  for  example,  a  ten-faced 
convex  shape,  we  select  seven  points  on  a  sphere  with 
a  unit  radius.  The  coordinate  values  of  these  seven 
points  can  be  given  by  (x,, yh  z,)  =  (sin  0,  cos  cp,,  sin 
0,  sin  cp,:,  cos  0,  ),  where  i  =  1-7  indicate  the  seven 
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Fig.  8.  Comparison  of  phase  function  and  degree  of  linear  polar¬ 
ization  for  randomly  oriented  six-faced  convex  aerosol  particles 
with  aspect  ratios  of  1,  2,  and  1/2. 


Fig.  9.  Comparison  of  the  phase  functions  for  oceanic  aerosol 
particles  of  various  shapes.  For  spherical  particles,  i.e.,  a/b  =  1, 
a  power-law  size  distribution  is  employed  to  smooth  out  the  reso¬ 
nant  fluctuations.  The  size  parameter  used  is  xmax  =  10. 


points  and  6,  and  cp,:  are  the  zenith  and  azimuthal 
angles  of  the  points,  respectively.  One  can  construct 
a  ten-faced  convex  shape  that  is  enveloped  by  a  unit 
sphere  by  properly  connecting  the  seven  points.  To 
obtain  a  convex  shape  with  a  desired  size  and  aspect 
ratio,  we  can  stretch  the  coordinate  values,  using  the 
relationship  (xi7  y,,  z,:)  =  (ax,,  ay, ,  bzt ),  where  a  and  b 
are  the  semiaxes  of  the  spheroid  that  envelops  the 
convex  particle.  We  obtain  an  oblate  particle  if  a  > 
b  and  a  prolate  particle  if  a  <  b. 

To  construct  a  concave  particle  we  select  two 
right-triangular  pyramids.  The  base  of  one  of  the 
pyramids  faces  up,  and  that  of  the  other  pyramid 
faces  down.  In  addition,  the  distance  from  the 
apexes  to  the  centers  of  the  pyramids  is  selected  to 
be  unity.  The  two  pyramids  share  the  same  origin. 
Thus  the  apexes  of  the  two  pyramids  are  confined 
on  a  unit  sphere.  Employing  the  preceding  coor¬ 
dinate  stretching  procedure,  we  can  obtain  various 
aspect  ratios  and  sizes  for  the  concave  geometry. 

We  select  two  refractive  indices,  m  =  1.38  +  i3.9  X 
10~9  and  1.75  +  iO.44,  which  represent  aerosols  with 
oceanic  sources  and  soots,  respectively.  Figure  8 
shows  the  phase  functions  and  the  degrees  of  linear 
polarization  computed  for  various  six-faced  convex 
aerosols  with  aspect  ratios  a/b  =  1/2,  1,  2.  The 
particles  are  assumed  to  be  randomly  oriented  in 
space.  For  a  size  parameter  of  5,  substantial  differ¬ 
ences  are  noted  for  the  three  aspect  ratios  in  the 
side-scattering  directions  about  120°.  For  a  size  pa¬ 
rameter  of  10,  the  differences  in  the  phase  function 


are  mainly  in  forward  scattering.  The  volume  asso¬ 
ciated  with  a/b  =  1  is  larger  than  those  for  a /b  =  1/2, 
2.  Thus  the  particles  with  a/b  =  1  have  a  larger 
scattering  capability  and  produce  a  stronger  forward 
peak.  However,  the  particle  volume  seems  not  to  be 
a  dominant  factor  in  determining  the  scattering  prop¬ 
erties.  For  example,  the  phase  functions  for  a/b  = 
1/2,  2  are  quite  different,  although  the  corresponding 
particle  volume  is  similar  in  these  two  cases.  The 
pattern  of  degrees  of  linear  polarization  is  not  sensi¬ 
tive  to  the  aspect  ratio,  except  for  a  size  parameter  of 
5  with  strong  absorption. 

Figure  9  shows  the  phase  functions  for  convex  and 
concave  aerosols  for  three  aspect  ratios  in  comparison 
with  those  for  spheroids  that  have  the  same  aspect 
ratios.  The  scattering  properties  of  spheroidal  par¬ 
ticles  were  computed  by  the  T-matrix  method.19  For 
a/b  =  1,  the  phase  functions  for  nonspherical  parti¬ 
cles  are  substantially  smaller  than  those  for  spheres 
in  the  scattering  angular  region  120°-180°;  in  partic¬ 
ular,  the  differences  in  the  results  for  spherical  and 
nonspherical  particles  are  pronounced  in  the  back- 
scattering.  For  a/b  =  2, 1/2,  we  also  see  substantial 
differences  between  smooth  spheroids  and  irregular 
convex  and  concave  particles. 

Finally,  we  compare  theoretical  phase  function 
results  with  experimental  data,  as  shown  in  Fig.  10. 
The  theoretical  results  correspond  to  an  ensemble 
of  randomly  oriented  concave  and  convex  particles. 
The  experimental  data  are  determined  from  the 
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Fig.  10.  Phase  functions  measured  by  the  microwave  analog  tech¬ 
nique  and  computed  from  the  FDTD  method  for  randomly  oriented 
convex  and  concave  particles  with  a  refractive  index  m  =  (1.5  + 
i0.005)  and  size  parameters  ranging  from  5.9  to  17.8. 


measurements  made  by  Zurell,20  who  used  the  mi¬ 
crowave  analog  technique  for  randomly  oriented 
convex  and  concave  particles  with  a  refractive  index 
m  =  1.5  +  i 0.005  and  size  parameters  ranging  from 
5.9  to  17.8.  The  data  have  been  reanalyzed  by  Pol¬ 
lack  and  Cuzzi.21  The  equivalent  spherical  results 
are  based  on  the  volume  of  the  particle  geometry 
used  in  the  theoretical  computation.  The  phase 
function  for  nonspherical  particles  computed  for  the 
FDTD  method  appears  to  match  the  experimental 
results. 

5.  Conclusions 

We  have  improved  the  Cartesian  FDTD  scheme  by 
properly  evaluating  the  dielectric  constant  for  par¬ 
tially  empty  cells  and  the  cell  edges  located  at  the 
interface  of  free  space  and  the  particle  medium.  We 
have  also  increased  the  efficiency  of  the  FDTD  pro¬ 
gram  by  employing  the  PML  technique  and  coding  in 
fortran90. 

Based  on  the  improved  FDTD  program,  the  accu¬ 
racy  of  various  approaches  to  dealing  with  the  sub¬ 
grid  variation  and  discontinuity  of  a  dielectric 
constant  at  the  particle  surface  was  investigated. 
The  magnitude  of  the  staircasing  approximation  of 
particle  shape  on  a  Cartesian  grid  is  inversely  pro¬ 
portional  to  the  ratio  of  the  grid  size  to  the  particle 
dimension.  For  a  small  size  parameter,  the  stair¬ 
casing  effect  is  substantially  large.  We  have  illus¬ 
trated  that  the  inverted  Maxwell-Garnett  rule  is  the 
most  effective  in  reducing  the  staircasing  effect. 
With  an  increase  in  the  size  parameter,  the  staircas¬ 
ing  effect  decreases,  and  the  method  that  one  uses  to 
evaluate  the  mean  permittivity  becomes  less  impor¬ 
tant.  We  also  investigated  the  effect  of  the  discon¬ 


tinuity  of  permittivity  at  the  interface  of  free  space 
and  the  particle  medium  that  exists  regardless  of  the 
pseudoparticle  shape  defined  in  the  Cartesian  grid. 
This  effect  must  be  properly  accounted  for  because 
the  electromagnetic  boundary  condition  implied  in 
the  computation  of  the  near  field  essentially  depends 
on  the  permittivity  at  the  surface  cell  edges.  Nu¬ 
merical  results  have  revealed  that  the  effective  per¬ 
mittivity  given  by  the  average  of  the  permittivity 
values  of  four  adjacent  cells  can  substantially  reduce 
the  discontinuity  effect. 

We  applied  the  improved  FDTD  code  to  the  study  of 
light  scattering  by  nonspherical  aerosols.  We  con¬ 
structed  convex  and  concave  particle  shapes  to  re¬ 
semble  dust  and  some  irregular  aerosols  and  carried 
out  calculations  of  phase  function  and  linear  polar¬ 
ization  patterns.  Comparisons  with  results  com¬ 
puted  from  the  T  matrix  for  spheroids  and  the  Mie 
theory  for  spheres  were  made.  We  showed  that  sub¬ 
stantial  differences  occur,  particularly  in  backscat- 
tering  directions.  Finally,  we  found  that  the 
theoretical  phase  function  computed  for  the  convex 
and  concave  particles  is  in  reasonably  good  agree¬ 
ment  with  available  microwave  analog  experimental 
results. 
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